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Using Local Operator Fluctuations to Identify Wave Function Improvements 

Kiel T. Williams and Lucas K. WagneiQ 
Department of Physics; University of Illinois at Urbana-Champaign, Urbana, IL 

A method is developed that allows analysis of quantum Monte Carlo simulations to identify errors 
in trial wave functions. The purpose of this method is to allow for the systematic improvement of 
variational wave functions by identifying degrees of freedom that are not well-described by an initial 
trial state. We provide proof of concept implementations of this method by identifying the need for a 
Jastrow correlation factor, and implementing a selected multi-determinant wave function algorithm 
for small dimers that systematically decreases the variational energy. Selection of the two-particle 
excitations is done using quantum Monte Carlo within the presence of a Jastrow correlation factor, 
and without the need to explicitly construct the determinants. We also show how this technique can 
be used to design compact wave functions for transition metal systems. This method may provide a 
route to analyze and systematically improve descriptions of complex quantum systems in a scalable 
way. 


I. Introduction 


First principles quantum Monte Carlo calculations [ij 
for solids are a promising way to go beyond density 
functional theory (DFT). These methods directly sim¬ 
ulate electron-electron correlations and can obtain very 
high accuracy on challenging materials using cur¬ 
rent state of the art techniques like fixed node diffu¬ 
sion Monte Carlo (DMC). Despite this success, the DMC 
method’s accuracy is limited by the fixed node approxi¬ 
mation, which allows for polynomial scaling of the com¬ 
putational cost with system size, but results in a DMC 
energy that is only an upper bound to the true ground 
state energy. In practical calculations, improvement of 
the accuracy and efficiency of fixed node diffusion Monte 
Carlo is reliant on improving trial wave functions which 
determine the fixed nodal surface. 

In order for a trial wave function to be appropri¬ 
ate for quantum Monte Carlo calculations, it should 
be compact and efficient to calculate. For application 
to bulk materials, it must also be size-extensive; that 
is, the total energy must scale with the system size. 
By far the most common trial wave function is the 
Slater-Jastrow wave function[^, Q, which is simple, ex¬ 
tensive, and initial guesses are easily obtainable from 
DFT codes. While truncated determinant expansions 
can be effective in describing small molecules (8|, they 
cannot be used in bulk materials because they are not 
size-extensive. Backflow wave functions^ have proven 
effective in homogeneous [l^ and inhomogeneous [ll| sys¬ 
tems, but worsen the computational scaling with system 
size of QMC methods and sometimes do not offer im¬ 
provement of results. It is thus of great interest, given 
a Slater-Jastrow wave function, whether there is a com¬ 
pact wave function that describes the most important 
improvements relative to the ground state. 

In this article, we present some first steps towards 


‘Electronic address: lkwagner@illinois.edu 



FIG. 1: Visual representation of the path through Hilbert 
space from the initial to exact wave function taken by the ex¬ 
act (black) and mimicked (red) projection operators respec¬ 
tively. 

a method that uses fluctuations of the local energy 
iLtI'(R)/d'(R), not to optimize a given parameterization, 
but to identify directions in Hilbert space that can im¬ 
prove trial wave functions. We first provide a summary 
of the principle of using the imaginary time projector 
exp(—riJ) to improve wave functions and the notation 
that will be used in the article. Then we show a proof of 
concept for multi-Slater Jastrow wave functions, in which 
this method is used to select determinants in the wave 
function. Finally, we show how the local energy fluctu¬ 
ations can be used to determine a priori what terms to 
add to a variational wave function for a transition metal 
system TiO. These results set the stage for data min¬ 
ing of many-body wave functions to determine how they 
should be improved. 

II. Theory 

In this work, we use ideas that have been known 
for a long time for optimizing parameters in wave 
functions |12 - ll4l | . but we follow more the work of 
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Holtzmannfisjl in that we would like to use the Feynman- 
Kac formulae to discover which parameterization to add 
to a given wave function. The quantum variational prin¬ 
ciple states that for any appropriately normalized trial 
wave function P), where R is the many-body elec¬ 

tron coordinate and P is a set of parameter values, the 
expectation value of the Hamiltonian of the system in 
state equals or exceeds the ground state energy of 
the Hamiltonian: 


Et > Eg, 


( 1 ) 


This minimization procedure provides an estimate of the 
set of associated {at} operator amplitudes. We define 

the local operators Hfc(R) = and a local energy 

Pl(R) = By expanding the projection operator 

to first-order and minimizing the square deviation, we 
find that: 


dk — T 


P^r(R) ifc^(R) 
\I/t(R) 'I't(R) 


'I'y(R)dR. 


( 8 ) 


where: 


Et{P) = (2) 

We then minimize Ej-{P) with respect to the parameter 
set P. Once this is done, we must alter the parameteriza¬ 
tion of the trial wave function to obtain further improve¬ 
ment in the energy estimate. Our ultimate goal will be 
not to optimize the parameters within a fixed set P, but 
to identify new parameters that must be added to P to 
improve the qualitative structure of the particular trial 
state. 

Iteratively applying the projection operator to a trial 
function 'I't’ produces a sequence of new wave functions: 

\^HiT))=e-^^\^T). (3) 

This converges to the exact ground-state wave function 
I^Gs) ia the infinite limit: 

lim |^'i/(r)) = I^'gs), (4) 

T—>-00 

Performing this operation directly corresponds to a pro¬ 
jector Monte Carlo method, such as diffusion Monte 
Carlo. The challenge in doing this is that compact rep¬ 
resentations of the operator exp(—riJ) are generally not 
known, and so the imaginary time dynamics must oper¬ 
ate in very high dimensions. Our objective here will be to 
find a compact representation of the short-time projector 
operator. 

We begin by considering an arbitrary set of linear op¬ 
erators {Ai}. Iteratively applying this set of operators 
generates a new sequence of wave functions: 


I^a) 


1 + y ^ ttiAi 


I'I't)- 


(5) 


For brevity, we define: 

\^aA=M'^t)- (6) 

We force the minimal set of operators {Ai\ to mimic the 
projection operator by minimizing the square deviation 
of A from T//: 


^ = -((F;i(R)-(R))i.(R)), (9) 

where we have assumed that elements of the set 
are orthogonal, or approximately orthogonal, such that 
Ri 5ik- Fig. [T] depicts this scheme pictori- 
ally, with the exact and mimicked projection operators 
represented by the black and tangential red arrows re¬ 
spectively. We see then that the mimicked projection 
operator evaluated for r = 0 can be viewed as a lin¬ 
earized approximation to the exact path to the ground 
state through Hilbert space. In this way, our approxima¬ 
tion to the projection operator identifies the most signifi¬ 
cant elements of Hilbert space absent from an initial trial 
state. 

The derivation of our method is similar in spirit to 
the stochastic reconfiguration (SR) of Sorella [H, [H- 
[H. The energy fluctutation potential method (EFP) 
also shares some similarities with our technique in its 
focus on the correlation between the local behavior of 
the energy and some chosen operator B [IS El. A 
set of operators Ai is a good set if only a few terms in 
Eqn[2]are non-zero, while a set with many small values in 
Eqn [3 is not an efficient descriptor of the wave function 
improvement. 


III. QMC Methodology 

We first compute the single-particle Hartree-Fock (HF) 
orbitals for a molecular system. We obtain all orbitals 
using the GAMESS computational package BjUSI. Core 
electrons were replaced by the corresponding Burkatzki- 
Filippi-Dolg pseudopotential [2^ with triple-^ basis sets. 

We perform variational Monte Carlo with the QWalk 
computational package [IS- We begin with a trial wave 
function of the Slater-Jastrow form: 

T = exp(C/)Det[(;ii(rj)], (10) 

We use the linear method of Umrigar et. al. B; 113 
to optimize the Jastrow U. The form of the Jastrow 
correlation factor [/ is a function of the electron and ionic 
coordinates: 


U = y^u{rii.rji,rij), 

ijl 


(^-A(R)-^'ff(R))'dR. 


(7) 


( 11 ) 
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where i and j indices represent electronic coordinates and 
I represents ionic coordinates. The functions u are given 
by: 


u{rii,rji,rij) = afc(r,/)+ 

k 

^kTmiakirii)aiirji)+ (12) 

m klm 

ak{rji)ai(rii))bk{rij), 


where the ak and bk functions have the general form: 


ak{r) = - 


1 - zirjr^ut) 


■ Pz{rlrcut)'' 


(13) 


and z{x) is a polynomial chosen to smoothly go to zero at 
T = Tcut [13 ■ This form of the Jastrow factor explicitly 
incorporates three-body interactions between two elec¬ 
trons and an ion. 


where the normalization factor is given by: 


iV, = 



(18) 


The two particle operators in Eqn[T3]are used to evaluate 
Eqn[9]and generate a list of important determinants miss¬ 
ing from the initial wave function. Hence, we can select 
the determinants most important to the exact ground 
state without the need to first evaluate those determi¬ 
nants. The entire process of wave function generation is 
summarized as such: . 


1. Obtain single-particle orbitals from a HE calcula¬ 
tion. 


2. Optimize single-determinant Slater-Jastrow: 

(R|d/)=e^(’'>^')Det[<^,(r,)]. (19) 


IV. Determinant selection 


The set of second-quantized two-body cre¬ 
ation/destruction operators, 


■Aij^ki — cjfecjjC.i'iCij, (14) 

offers one possible choice of linear operators Ak in Eqn[51 
If i,j are occupied orbitals and k,l are unoccupied or¬ 
bitals, then applying a Aij^ki to a Slater determinant gen¬ 
erates an additional excited-state determinant. The ele¬ 
ments of the two-body reduced density matrix (2-RDM) 
are given by the expectation values of these two-body 
creation/destruction operators. We thus make the anal¬ 
ogy with local energy to define a local density matrix 
element, given a wave function Id/x): 


Pijkl (I^) 


Or, explicitly: 


Aij^ki'^T(R-) 

'I'ji (R) 


(15) 


Pijki(R)=J2 f 4'ki'r'aWii'rb) 

(16) 

X (/)i (ra) (rfe) 4 - 5 . (i?"f,) [R)dr'^ dr[,. 

where R = (ri, r 2 ,..., tat), = 

(ri, r 2 ,..., r(j, ..., r^r) refers to the set of coor¬ 
dinates generated by changing the positions of two 
electrons, and we have omitted overall normalization. 
We evaluate this 2-body integral in a QMC calculation 
by sampling the coordinates and r(, from the sum over 
orbitals /(r) = X)i?^?(''') the many-body electron 
coordinate R from 'I'^(R) [i^. With this, the expression 
given in Eqn[Tn]can be rearranged to give: 


Pijkl (I^) — 


1 

N.NjNkNi 



'i'(R) 


(l>kK)<l^iiK)M'ra)(l)jirb)' 


3. Rank 2-RDM elements by covariance of 

4. Add corresponding determinant to the expansion: 

|4'new) = |4'oid)+aiet^[c4c|;Cf-iC^j]|Hartree-Eock). (20) 

5. Optimize coefficients {ci} of [d^new) using the linear 
method. 

Iterating steps 3-5 of this process generates a determi- 
nantal expansion of arbitrary length up to the full size of 
the active space. 


A. H2 molecule 

For the case of H 2 , we restrict our active Hilbert space 
to the set of bonding/antibonding tr-symmetry orbitals. 
Fig. shows the contours of the sampled amplitude 
(^) the local operator associated with a 2- 
body b ^ a bonding-to-antibonding excitation in an iso¬ 
lated hydrogen dimer versus the sampled local energy 
ElCR.) for each of two trial states: 

'^SJ = e^<))bt(D)<)>bi(c2) 

'^MSJ = {ci(j)bj{ri)(j)bl -f C2()<at(D)<)>a4,), 

where ^sj and 'I'msj are the Slater-Jastrow and multi- 
Slater Jastrow wave functions containing the bonding (j)b 
and antibonding (pa single-particle orbitals respectively. 

The line segments on each panel in Fig [2] indicate the 
principal components of the resulting distribution. These 
components are given by the eigenvectors of the covari¬ 
ance matrix of the local energy distribution taken with 
respect to the local operator pa.j.aj^{,.|.{,j^(R): 

f ^p,p ^p,El \ 

\<^El,p <^El,ElJ 


fK)f{r'b) 


( 17 ) 
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El (R) (Hartree) El (R) (Hartree) 

FIG. 2: Amplitude of 2-body b ^ a bonding-to-antibonding 
excitation 6 .^ 64 ^ (R.) versus local energy El(R.) for two 

different trial wave functions, with corresponding principal 
components of the distribution indicated. The Slater-Jastrow 
wave function used to generate the left panel did not include 
the CSF corresponding to this 2-body excitation, while the 
wave function used to generate the right panel does. The 
principal components rotate upon the addition of this CSF. 


in this two-dimensional representation. The rotation of 
the principle components relative to the axes in the left 
panel of Fig. [5] shows that the covariance matrix con¬ 
tains nonvanishing off-diagonal elements. It follows that 

are correlated for this single¬ 
determinant trial state. After the addition of the asso¬ 
ciated 6 —a determinant to the wave function in the 
right panel of Fig. [2l the principle components rotate to 
align with the axes, indicating that the covariance ma¬ 
trix has become diagonal. This implies that the covari¬ 
ance between the local energy El(R) and local opera¬ 
tor Pa.fa 4 .f)i.b 4 . (R) has vanished, and the two variables now 
have zero covariance. That is, a key element absent from 
the initial trial state has been identified and added based 
on the covariance, pushing the wave function closer to 
the exact ground state. 


B. Dimer Molecules 

As a further proof of concept, we apply the covari¬ 
ance method to select determinants for a set of stretched 
molecules: H 2 (0.88 A bond length), N 2 (1.7 A bond 
length), O 2 (1.6 A bond length), and F 2 (1.5 A bond 
length). By stretching the molecules, the electron cor¬ 
relations are enhanced, increasing the strength of the 
covariance signal. We obtain single-particle orbitals 
for each system from a restricted open-shell Hartree- 
Fock (ROHE) calculation using GAMESS. This method 
doubly-fills molecular orbitals (MOs) to the greatest ex¬ 
tent possible, and places remaining unpaired electrons 
into singly-filled MOs. We limit our active space to a set 
of bonding and antibonding MOs with cylindrical sym¬ 
metry and either a- or 7r-symmetry. Other states exist 


within the full orbital space, but their inclusion yields 
only small improvement to the final wave function and 
system energy. Because different methods of determinant 
selection produce significantly different rates of energy 
convergence [s^, the covariance-based method we have 
described can yield interesting results even at the level 
of a multi-Slater-Jastrow ansatz. Our chief objective in 
this section is to show that the covariance technique can 
select the most significant determinants for a particular 
molecule before performing a variational optimization of 
the wave function. 

We consider only 2-particle excitations featuring 1 par¬ 
ticle in each spin channel. We compare these results to 
those obtained with the usual configuration interaction 
method with singles and doubles excitation (CISD). This 
is natural for molecules such as N 2 with a ground state 
singlet spin configuration, though it can lead to the exclu¬ 
sion of significant excitations in molecules like O 2 which 
contain unpaired electrons. Fig. [3] compares the nor¬ 
malized weight of each CSF in conventional CISD, the 
optimized weight of each CSF in a multi-Slater-Jastrow 
ansatz, and the local energy covariance for each relevant 
CSF in each material respectively. We see that the deter¬ 
minant orderings predicted by both traditional CISD and 
our method based on local energy covariance are equiv¬ 
alent for each system across the dominant particle exci¬ 
tations. This indicates that the path to the ground state 
through Hilbert space obtained by successively applying 
the projection operator is approximately equivalent to 
that produced by the usual Cl procedure in this case. 

From Eqn[Hl we see that the covariance signal in a 2- 
RDM element should fall identically to zero once the cor¬ 
responding excitation has been added to the trial state. 
In practice, we observe that the signal in an added ex¬ 
citation falls significantly once it has been added to the 
trial wave function, but it does not vanish entirely. This 
is a consequence of the Jastrow factor in the trial state, 
which spoils the orthogonality of determinants and pre¬ 
vents the covariance in each 2-RDM element from fully 
vanishing. That is, the overlap matrix Stk is non-diagonal 
when a Jastrow factor multiplies the determinant expan¬ 
sion. Practically speaking, this slight non-orthogonality 
did not seem to affect the performance of the technique. 
Assuming that Sik is diagonal allows the method to scale 
linearly with the number of excitation operators, so we 
use this approximation. 

We also find the rate of energy convergence for the 
predicted CSF ordering in each model molecular system. 
Fig. [4] shows the variational Monte Carlo energy of an 
optimized multi-Slater-Jastrow wave function as a func¬ 
tion of the CSFs included in the trial state. The CSFs 
are ordered here according to the weight given by a con¬ 
ventional CISD calculation. We see that the energy con¬ 
verges rapidly with respect to the number of CSFs in¬ 
cluded in the wave function. This explicitly illustrates 
that the CISD method and our covariance-based tech¬ 
nique can drive the initial trial state asymptotically close 
to the exact ground state. 
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FIG. 3 : Comparison of normalized signal strength for dif¬ 
ferent estimators of relative CSF importance for stretched 
dimers of H2, N2, O2, and F2 respectively. The CSFs are 
arranged such that the optimized final CSF weight declines 
monotonically from top to bottom. Each indicated excitation 
is a 1- or 2-particle excitation that includes both itself and any 
symmetry-related partners. For example, (ni,Tri) {Tii,'Jn) is 
a 2-particle excitation that excites a bonding yr-orbital elec¬ 
tron to an antibonding tt* orbital of the same angular mo¬ 
mentum [x or y) in each spin channel. On the other hand, 
{TTi,TTj) —>■ (TTjjTTi) involves a two-body exchange. 


FIG. 4 : Added spin-up/spin-down CSF excitations vs. associ¬ 
ated variational Monte Carlo energy in a multi-Slater-Jastrow 
wave function for the CSF ordering suggested by conventional 
CISD for each considered model system. 


Finally, we also assess the degree to which the covari¬ 
ance in a 2-RDM element predicts the energy gain ob¬ 
tained from adding the associated determinant to the 
trial state. Fig. [5] compares the decrease in total system 
energy obtained from each additional CSF with the cor¬ 
responding covariance signal. We observe that the energy 
gain and the covariance signal are negatively correlated 
with one another. This correlation indicates that the co- 
variance in a 2-RDM element can be used as a proxy for 
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FIG. 5: The normalized covariance signal of CSFs versus the 
decrease in energy obtained from adding a CSF to the trial 
state. Significant negative correlations exist between the two 
values. The shading is provided as a visual guide. 


estimating the energy change from adding a determinant 
to the trial state. 

As a method of determinant selection for these sys¬ 
tems, this technique is less efficient than using Cl to 
determine the weights, and the results are similar. We 
therefore would not recommend this technique as a se¬ 
lection method for small molecules. However, the point 
of this section is that the energy fluctuations can be data 
mined to find the correct directions in Hilbert space to 
improve trial wave functions. In the case of stretched 
dimers, it is well-known that the most important im¬ 
provement over Slater-Jastrow consists of multiple de¬ 
terminants, and the energy fluctuation technique selects 
the correct ones. 


V. Comparing real and orbital spaces: TiO 
molecule 
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FIG. 6: (a.) The covariance of the 1-RDM with the local 
energy El for TiO in the J12 wave function, (b.) The co- 
variance of the pair distribution ^(r) with the local energy 
El for the J12, J30, and JS30 wave functions in both spin 
channels (right). J12: Slater-Jastrow state with 12 parame¬ 
ters per atom in the 3-body part of the Jastrow factor; J30: 
J12, but with 30 3-body terms instead of 12; JS30: J30, but 
with a spin-dependent 2-body portion of the Jastrow factor. 



J12 

JS12 

J30 
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Number of CSFs 


We now proceed to use the technique to selectively 
improve wave function parameterizations in a more chal¬ 
lenging case. As an example of a system where we do 
not know a priori the most important degrees of free¬ 
dom, we consider a transition metal molecule, TiO. The 
dynamic correlation present in transition metal systems 
is larger than in s-p systems like the dimers considered 
above, so the Jastrow factor could be expected to play a 
larger role mill. 

In this case, the 1-RDM can be used to select the or¬ 
bital space by calculating the covariance of the diagonal 
elements of the 1-RDM with the local energy. In the case 
of the stretched N 2 dimer, the elements corresponding to 
the bonding and antibonding orbitals have a covariance 
with the local energy of approximately 0.001, while other 
orbitals have much smaller signals. 

In Fig [HI the covariances of the 1-RDM and the real- 


FIG. 7: Total TiO VMG energy vs. the number of included 
GSFs using different Jastrow factors. Note that the decline 
in energy is quite modest with respect to the number of in¬ 
cluded CSFs, but falls dramatically when spin-dependence is 
incorporated in the 2-body portion of the Jastrow factor. 


space electron-electron correlation function g{r) distance 
are shown. The covariance signal for the 1-RDM is very 
small, and in fact the selection of determinants was very 
difficult in this case, although we do obtain larger sig¬ 
nals for the p and d states as one would expect. On 
the other hand, for our starting wave function, labeled 
J12, with 12 three-body parameters per atom, there is 
a large spin-dependent covariance with g{r). Increasing 
the number of parameters in the Jastrow factor with- 
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out making it spin-dependent (J30, with 30 3-body pa¬ 
rameters per atom) does not resolve this covariance, but 
adding in four spin dependent parameters (JS30) imme¬ 
diately reduces the covariance signal to nearly zero across 
all separation distances. 

Since the determinant selection of TiO via energy co- 
variance was not efficient, we used a Cl calculation with 
sextuple excitations into 8 virtual states to select CSFs, 
then formed a set of multi-Slater Jastrow wave functions. 
If the covariance analysis was correct, then we would ex¬ 
pect the spin-dependent terms in the Jastrow to be most 
effective in lowering the energy, followed by either the ex¬ 
tra three-body terms or multiple determinants. As can 
be seen in Fig 0 this supposition is correct: with only 
four parameters, the spin-dependent terms lower the en¬ 
ergy by nearly 10 mHartree, while 30 determinants or 
a similar number of 3-body parameters are necessary to 
achieve that decrease in energy. 

This example illustrates some the strengths and weak¬ 
nesses of this covariance-based selection. If the set {Ai} 
is selected in a basis that does not describe the needed im¬ 
provement efficiently, in this case the determinant basis, 
then it is not the best tool. On the other hand, if sev¬ 
eral different basis sets are used, then the best basis can 
be used to improve the wave function. In this case, we 
learned that a spin-dependent Jastrow factor can improve 
the energy significantly for magnetic molecules, while the 
determinant basis is not an efficient way to improve the 
wave function for this molecule. The cost for performing 
these calculations was about a factor of two larger than 
a variational Monte Carlo calculation and much smaller 
than the energy optimization technique. 

VI. Conclusion 

We have presented an outline of a technique to select, 
not just terms in a many-body ansatz, but which type of 


ansatz with which to proceed. For example, the selection 
method can quickly determined whether a determinant- 
type basis is appropriate by evaluating the 1-RDM co- 
variance with the local energy. Similarly, if an explicitly 
correlated approach such as a Jastrow is more appro¬ 
priate, then the covariance of the local energy with the 
electron-electron distance g{r) is large. The computa¬ 
tional cost of this assessment is quite low: g(r) is essen¬ 
tially zero cost over a VMC energy evaluation, and the 
1-RDM is approximately a factor of two additional, re¬ 
gardless of system size. This is much less expensive than 
attempting energy minimization on multiple ansatz. 

As proof of concept, we demonstrated that the se¬ 
lection technique both selects the correct directions in 
Hilbert within a defined ansatz space, and also can select 
between alternate viewpoints of the electron correlation 
problem. We demonstrated the former by selecting de¬ 
terminants for stretched dimer molecules, and the latter 
by differentiating between short range ’dynamic’ correla¬ 
tion best described by a Jastrow factor and long range 
’static’ correlation best described by multiple determi¬ 
nants in the transition metal oxygen system TiO. Using 
standard wave functions for this problem, the dynamic 
correlation in TiO is more important. This work forms 
the base for an algorithm in which the local energy can 
be analyzed directly in the many-body space using fea¬ 
ture extraction techniques to describe the most efficient 
basis in which to improve many-body wave functions. 

This material is based upon work supported by NSF 
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dation Graduate Research Fellowship Program under 
Grant Number DGE-1144245 (K.T.W.). We also ac¬ 
knowledge computer resources from the Campus Clus¬ 
ter program at Illinois. Useful conversations with David 
Ceperley are gratefully acknowledged. 
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